setwd("/home/uec-00/yapingli/code/Allele_specific_methylation/")
for (e in commandArgs(TRUE)) {
  ta = strsplit(e,"=",fixed=TRUE)
  if(! is.na(ta[[1]][2])) {
		if(ta[[1]][1] == "k"){
			k<-ta[[1]][2]
			
		}
		if(ta[[1]][1] == "s"){
			s<-ta[[1]][2]
			}
	}
}
		fileName=paste("methylCGsRich_ASM_AllSnp_",k,"Merge_chr",s,"_table.txt",sep="")
		x<-read.table(fileName,header=F,sep="\t")
		ASM<-array()
		i=1
		while(i<=length(x[,2])){
			temp<-x[i,1]
			y<-x[x[,1]==temp,]
			len<-length(y[,2])
			
			minRange<-min(y[,2])
			maxRange<-max(y[,2])
			region=paste(minRange,maxRange,sep="~")
			
			A<-c(sum(y[,3]),sum(y[,4]))
			B<-c(sum(y[,5]),sum(y[,6]))
			
			p<-fisher.test(rbind(A,B),workspace = 2000000)
			logP<-0-log10(p$p.value)
			methy_A<-sum(y[,3])/(sum(y[,3])+sum(y[,4]))
			methy_B<-sum(y[,5])/(sum(y[,5])+sum(y[,6]))
			
			methy_ratio=0
			ifelse(methy_B != 0,methy_ratio<-log2(methy_A/methy_B),methy_ratio<-c("NA"))

			thisRow<-cbind(temp,region,length(y[,2]),methy_ratio,p$p.value,logP)
			ASM<-rbind(ASM,thisRow)

			i<-i+len
		}

		pAdjust<-p.adjust(as.numeric(ASM[,5]),method="BH")
		pAdjustLog<-0-log10(pAdjust)
		ASM<-cbind(ASM,pAdjust,pAdjustLog)
		outputName=paste("methylCGsRich_ASM_AllSnp_",k,"Merge_chr",s,"_pValue_ASMblock.txt",sep="")
		write.table(ASM, outputName,sep="\t",quote=F,col.names=F,row.names=F)
	


